Strong Collapse Turbulence in Quintic Nonlinear Schrodinger Equation 
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We consider the quintic one dimensional nonlinear Schrodinger equation with forcing and both 
linear and nonlinear dissipation. Quintic nonlinearity results in multiple collapse events randomly 
distributed in space and time forming forced turbulence. Without dissipation each of these collapses 
produces finite time singularity but dissipative terms prevents actual formation of singularity. In 
statistical steady state of the developed turbulence the spatial correlation function has a universal 
form with the correlation length determined by the modulational instability scale. The amplitude 
fluctuations at that scale are nearly-Gaussian while the large amplitude tail of probability density 
function (PDF) is strongly non-Gaussian with power-like behavior. The small amplitude nearly- 
Gaussian fluctuations seed formation of large collapse events. The universal spatio-temporal form 
of these events together with the PDF for their maximum amplitudes define the power-like tail of 
PDF for large amplitude fluctuations, i.e., the intermittency of strong turbulence. 

PACS numbers: 47.27.-1, 42.65.Jx, 52.38.Hb 



I. INTRODUCTION 

A nonlinear Schrodinger equation (NLS) 

zVt + VV + a|V'|V + /3|V'lV = 0, (1) 

describes a wide class of interacting nonlinear waves and 
Bose-Einstein condensates. Here t is the time, the Lapla- 
cian is considered in the general dimension D, the 
constants a and /? correspond to the cubic and quintic 
nonlinearities, respectively. Generally, ajV'P'i/' giives the 
leading order nonlinear interaction. In many nonlinear 
systems a can vanish, which results in the quintic non- 
linear Schrodinger equation (QNLS) 



#t-HVV + /3|^lV = 0. 



(2) 



QNLS occurs e.g., in the Bose-Einstein condensate where 
the s-wave scatteringlength is set to zero by tuning Fes- 
hbach resonance [H, H- QNLS also occurs for general 
NLS type system near the transition from supercritical 
to subcritical bifurcations 0, 3], pattern formation (in 
the context of quintic Ginzburg-Landau equation if /3 is 
the complex constant) ^ and dissipative solitons (e.g., 
in lasers) [6]. Another possible experimental realization 
of equation ^ is the optical pulse propagation in optical 
fiber using a nonlinear compensator of nonlinearity Q- 

The standard cubic NLS (equation ([1]) with /3 = 0) is 
integrable in dimension one (ID) by the inverse scattering 
transform f8!| with global existence of all solutions. In 
contrast, QNLS ([2]) with positive real /3 and any D > 1 
can develop a finite time singularity (blow up) such that 
the amplitude of solution reaches infinity in a finite time. 
The blow up is accompanied by dramatic contraction of 
the function ip spatial extent, that is called wave collapse 
or simply collapse 0, ■ A sufficient condition for the 
collapse is 77 < 0, where 



is the Hamiltonian (energy) and the equation ^ can be 
rewritten in the Hamiltonian form 



SH 



(4) 



The case D = 1, which we consider below, is critical 
because any decrease of the power of nonlinearity in ^ 
(i.e., replacement of by 7 < 0) results in 

the global existence of the solutions |lll - [l3| for any real 

Collapse of QNLS is not physical and near singularity 
different physical regularization mechanisms come into 
play. These can be numerous nonlinear dissipation mech- 
anisms such as inelastic collisions in Bose-Einstein con- 
densate which results in loss of particles from the conden- 
sate [T|, optical breakdown and formation of plasma in 
nonlinear optics [l^ . or numerous non-dissipative regu- 
larization effects such as nonlinear saturation in laser- 
plasma interactions [l^, different dispersive effects or 
non-paraxiality of optical beam (see e.g., [l6j). 

In this paper we consider ID QNLS with linear and 
nonlinear dissipation so that QNLS ^ is replaced by 
the following regularized QNLS (RQNLS): 



i-0t + (1 — iae)dl'>p + (1 + zce)|-0|'*V' = ie0, 



(5) 



(3) 



which can be also called as a complex quintic Ginzburg- 
Landau equation. Here x is the spatial coordinate re- 
placing general r, < e ^ 1 is a small parameter so that 
to the leading approximation e — 0, QNLS is valid. 
The coefficient a ~ 1 determines linear dissipation and 
the coefficient c ~ 1 is responsible for nonlinear dissipa- 
tion. The linear dissipation has a viscosity-like form and 
can be resulted, e.g., from angular-dependent losses or 
the optical filtering [13, [13 ■ The nonlinear dissipation 
in RQNLS corresponds to the three-photon absorption 
in optics [l3| or four-body collisions which cause loss of 
atoms from the Bose-Einstein condensate UM ■ The term 
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ie(j) describes the general forcing in the system. The spe- 
cific examples of the realization of RQNLS ([5]) are e.g., 
a propagation o f lig ht in a ring cavity with Kerr nonlin- 
earity (see e.g., [20]) or any quite general propagation of 
waves in nonlinear media with complex dispersion and 
nonlinear dissipation pll |. 

The right hand side (rhs) of equation ([S]) provides forc- 
ing and depends on the specific physical model. We con- 
sider two types of forcing. First is a deterministic forcing 

</- = (6) 

which corresponds to a linear instability (amplification) 
in a system. Here b is the linear integral operator over 
X such that its spatial Fourier transform bk is the mul- 
tiplication operator (f>k = bj-ipk, where bk determines 
/c-dependence of the amplification. Below, if not men- 
tioned otherwise, we implicitly assume the simplest case 
of A:-independent amplification as bip = bip and thus, the 
equation ^ takes the following form 

<^ = 6^, (7) 

where b is the positive constant. 

The second type is a random additive forcing 

(t>^at,x), (8) 

where zero in average stochastic term ^ is a random Gaus- 
sian variable which is ^-correlated in time and has a finite 
spatial correlation such that 

{atl,Xi)C{t2, X2)) = S{t, - h)x{\xi - X2\). (9) 

Here 

x{\xi - X2\) = 6gexp[-|xi - X2\/lc], (10) 

Ic is the correlation length of pump, (. . .) denotes aver- 
aging over the statistics of ^, and bg is the normalization 
constant. 

Forcing results in the pumping of energy into the sys- 
tem described by RQNLS and subsequent formation of 
multiple collapse events randomly distributed in space 
and time as shown in Figure[T^. Figure[TlD shows a zoom- 
in into a temporal evolution of a spatial profile of a typi- 
cal collapse event which involves growth and subsequent 
decay of collapse amplitude. 

After the initial transient, the solution of RQNLS 
achieves a statistical steady state (i.e., state of the de- 
veloped turbulence) as shown in Figure [5] through the 
time dependence of the following integral 

J iV'pdx, (11) 

which has a meaning either of the number of particles in 
Bose-Einstein condensate, or an optical power (or some- 
times energy) in optics as well as it is called by a wave 
action in oceanology and many other nonlinear wave ap- 
plications Below we refer to N as the number of 
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FIG. 1; (Color online) Spatial-temporal form of solution of 
RQNLS ® with the deterministic forcing Q, b = 10*, 
e — 2 X 10~^, and a — c — 1 for different moments of time, 
(a) Collapses (seen in plot as sharp peaks) occur at random 
spatial positions and random times, (b) Zoom-in at evolution 
of single collapse event is shown for the smaller time interval 
and the smaller spatial extent compare with (a). It is seen 
that the collapse amplitude grows initially, then goes through 
the maximum and finally decays due to dissipative effects. 
Colors are added to (a) and (b) to help distinguish different 
moments of time. 



particles. It is seen in Figure ^jp that the dissipation 
is important near large collapses while forcing works all 
time (because forcing is oc N). In the statistical steady 
state the pumping of particles (forcing) in average is com- 
pensated by the dissipation which insures the state of the 
developed turbulence. 

The focus of this paper is to describe the strong tur- 
bulence in RQNLS characterized by these nearly-singular 
collapse events. By strong turbulence (we also call it 
strong collapse turbulence) we mean the turbulence with 
strong non-Gaussian fiuctuations as opposed to the weak 
turbulence with nearly Gaussian fluctuations plj . The 
strong non-Gaussian fiuctuations are usually refer to as 
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FIG. 2: (Color online) Time dependence of the number of par- 
ticles A*' for numerical simulation of RQNLS (O with the de- 
terministic forcing (O (left axis, thick solid line) superimposed 
with the time evolution of maximum value of (right axis, 
solid line). Parameters are 6 = 10*, a = c= l,e = 2x 10~^ 
and the random small amplitude initial condition is used, (a) 
Initially A'^ grows until the statistical steady-state of devel- 
oped turbulence is achieved so that A'^ remains constant in 
average over time, (b) Zoom-in at the smaller scale in t. 
Solid thick curve (blue) corresponds to N{t). Solid thin curve 
(green) shows time dependence of maximal spatial value of 
\^p\. Scale for A'^ is on the left vertical axis and scale for max 
is on the right vertical axis, respectively. It is seen that sharp 
decreases of A' is due to the dissipation from large collapses. 
These decreases are compensated in average by growth of A'' 
from the linear forcing. 



intermittency of turbulence [2^. The classical example 
of the strong turbulence is the Navier-Stokes turbulence. 
The old idea of the description of strong turbulence in 
the Navier-Stokes equations through singularities of the 
Euler equations still remains unsolved 22]. The forced 
Burgers equation is a very rare example of an analytical 
description of strong turbulence in which the tail of the 
probabihty density function (PDF) for negative gradients 
follows a well-established (—7/2) power law [2^, domi- 
nated by the spatio-temporal dynamics near formation 
of singular shocks (i.e., by pre-shocks). The spectrum 
of the strong optical turbulence was considered for the 
two-dimensional cubic nonlinear Schrodinger equation in 
Ref. [l^l and power law tail of PDF of the amplitude 
fluctuations was considered in Refs. (2^ and [26|. Strong 
turbulence in RQNLS was first studied in Refs. [l^l and 
[2^ with the (—8) PDF power law scaling suggested in 
Ref. jl^l- Here we show that (—8) is only an approximate 
scaling and it is determined by fluctuations of the waves 
at background which seed collapses as well as by the self- 
similar form of collapsing solutions. We also found that 
the power law of PDF tail is only weakly sensitive to 
the type of forcing (linear ampliflcation ([T]) vs. additive 
random forcing ([8])) for e <C 1 showing the universal tur- 
bulent picture. 

The paper is organized as follows. In Section |ll] we 
consider the small amplitude fluctuations of the back- 
ground of RQNLS turbulence, and show the universality 
of the spatial and temporal correlation functions and re- 
late the correlation scale of these fluctuations to the scale 
of the modulational instability. In Section IIIII we review 
the collapsing self-similar solution of RQNLS for e = 
and establish the universality of the self-similar solution 
for e 7^ as a basic building block for large amplitude 
fluctuations. In Section IIV Al we provide the results of 
numerical calculations of PDF for the amplitude fluctu- 
ations and for the collapse maximums. In Section IIVBI 
we derive the analytical expression for the tail of PDF 
for amplitudes and compare with numerics. In Section 
V we discuss the details of the numerical methods used. 
In Section VI the main results of the paper and future 
directions are discussed. 



II. MODULATIONAL INSTABILITY AND 
FLUCTUATIONS OF THE BACKGROUND 



The forcing term in right hand side (rhs) of RQNLS 
^ pumps the number of particles N into the system un- 
til the statistical steady state (also can be called by a 
developed turbulence state) is reached as shown in Fig- 
ure [2^ for the particular example of the deterministic 
forcing ([7]). The system in developed turbulence state 
does not have memory of initial conditions because of 
the modulational instability. That instability was de- 
rived for the cubic NLS [2^ (see also e.g. [2l[) but 
it is straightforward to generalize it for RQNLS as fol- 
lows. RQNLS ([S]) with e = has a spatially uniform 
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solution i/j = -00 exp(i|'0o|^i). Linearization on the back- 
ground of that solution in the form ip = exp(j|'0o|^<) [?Ao + 
(5-0 exp(z/i + ifcx)] with \6ip/tpo\ <C 1 gives the following 
instability growth rate v for the wavenumber k: 

The instability occurs for |fc| < 2|0oP with the unstable 
branch i?e(i/(fc)) = i^(fc) > reaching maximum of i^(fc) 
for 

= 21^0^ (13) 

For e 7^ 0, the expressions and are still ap- 
proximately valid provided v{kmax) = 2|V'o|'' ^ I/'''/: 
where r/ is the typical time of forcing. For the deter- 
ministic forcing Tf — {eb)^^ while for the stochas- 
tic forcing (|8]), Tf can be estimated from the condition 
that dissipation in average is compensated by the forc- 
ing. Then Tf ^ min[(ecpQ)~^, (ea/cQ)^^], where po is 
the typical amplitude of the fluctuation of Itpl and ko 
is the typical wavevector which is estimated from (1131) as 
ko ^ kmax- Recalling now our assumption that a ^ 1 
and c ~ 1 we arrive to a simpler estimate r/ ~ (epg)'^, 
e.g., using the parameters of Figure [5] we obtain that 
l/T/~20<zy(fc„,a.)^3-102. 

Thus the dynamics of the background of turbulence 
can be characterized by the typical amplitude of the fluc- 
tuations Po and spatial scale 1/kmax = l/(2^/^Po)- For 
simulations we define po = (iV/Lo)^^^, where Lq = J dx 
is the computational domain (without loss of generality 
we set Lg = 1 in our simulations as described in Section 
V) . We determine a correlation length Xcorr of ip through 
the full width at half maximum (FWHM) of the spatial 
correlation function S{x) = (0(xo, i)V'*(a;o +x)). Here, 
(...) denotes averaging over space and time which is as- 
sumed by the ergodicity to give the same result as the 
average over the ensemble of simulations. Figure [3] shows 
\S{x)/ S{0)\ vs. x/xcorr for a set of simulations with dif- 
ferent sets of parameters. Each curve is calculated in 
the statistical steady state. It is seen that |S'(a;)/5'(0)| 
is well approximated by the universal function of x/Xcorr 
which is close to a Gaussian function exp [—(x/xq)'^] while 
5(0) = Pq and xq is chosen from the condition to have the 
same FWHM of the Gaussian function and |S'(a;)/5(0)|. 
It implies xo/xcorr — l/(2-\/ln2). Note with the increase 
of the ensemble of simulations (i.e., the simulation time 
used to calculate S{x)), the value of Im(5(x)) approaches 
to 0. 

We also determine a correlation time tcorr of ip through 
FWHM of the temporal correlation function P{t) = 
{'ip{x,to)ip*{x,to + t)). Figure H shows \P{t)/P{0)\ as 
a function of normalized time t/ tcorr calculated in the 
statistical steady state. It is seen that for different set 
of parameters, \P(t)/ P{Q)\ is well approximated by the 
universal function of t/ tcorr- The fluctuations in the tails 
of \P{t) / P{Q)\ in FigureUlare due to the finite size of the 
statistical ensemble. 
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FIG. 3: (Color online) Absolute value of the normalized 
spatial correlation function |S'(a::)/5'(0)| vs. x/xcorr- Pa- 
rameters are e = 5 • 10"'^, a = c = 1,6 = 10^ (solid 
blue line), e = 2 ■ 10"^, a = c = l,b = 10* (dashed- 
dotted blue line), e — 10~^,a = c — l,h — 10* (dashed 
blue line), e = 2 ■ 10"^, a = c = l,b = 2 • 10^ (dot- 
ted blue line), e = 2 • 10^^, a = l,c = 2, 6 = 10* (green 
crosses), and e = 2 • 10"^, a = 2,c = 1,6 = 10* (light 
blue circles). Red squares correspond to a Gaussian function 
exp [-41n2(a;/2;corr)^] which has a unit FWHM. 




FIG. 4: (Color online) Absolute value of the normalized 
temporal correlation function \P{t) / P{Q)\ vs. t/tcorr- Pa- 
rameters are e = 5 • 10"'^, a = c = l,fe = 10* (green cir- 
cles), e = 2 ■ 10"^, a = c = 1,& = 10* (dashed-dotted 
blue line), e = 10""*, a = c = l,b = 10* (dashed blue 
line), e = 2 ■ 10"^, a = c = 1,6 = 2 ■ 10^ (dotted blue 
line), e = 2 ■ 10~^,a = l,c = 2, 6 = 10* (red crosses), and 
e = 2 • 10~^ a = 2, c = 1, 6 = 10* (solid blue line). 
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FIG. 5: (Color online) Dependence of p§ on x'Jj.,., where Xcorr 
is the FWHM of spatial correlation function S{x). Circles 
correspond to the case of fc-independent deterministic forcing 
with a=l,b= 10^c = 1, e = 2 • 10"^ a = 1, 6 = 2 ■ 10^ c = 
l,e = 2 - 10-3; a = l,fe = 3 - 10^c = l,e = 2 ■ 10"^ a = 2,h = 
10'',c = l,e = 2 • 10"^ a = l,fe = 10^, c = 2, e = 2 ■ 10"^ 
a = l,b = 10^, c = 1, e = 5 • IQ-^, from leftmost to rightmost. 
Squares are for the "fc-limited" deterministic forcing ((6)1 , (|15|) 
with Kutoff = IOtt for a = 1,6 = 10^,0 = l,e = 2 ■ 10"^; 
a = 1,6 = 3 - 10^c = l,e = 2 ■ 10"^ a = 2, 6 = 10'',c = l,e = 
2 ■ 10^3; a = 1,6 = 10'',c = 2, e = 2 ■ 10"^ a = 1,6 = 10**, c = 
1, e = 2 ■ IQ-'', from leftmost to rightmost. Diamonds are for 
the random forcing with a = 1, 6g = 250, c = 1, e = 5 • 10"^; 
a = l,6g = 1250, c = l,e = 5 • 10"^; a = 1,6^ = 1250, c = 
l,e = 2-10-^ a = 1,63 = 2500,c = l,e = 2-10"^ a = 2,69 = 
5000, c = 1, e = 2 ■ 10-=^;a = 1, 6g = 5000, c = 2, e = 2 • 10'^ 
a = l,6g = 5000, c = l,e = 2 ■ 10"^, and Ic = 0.02 in all 
cases from leftmost to rightmost. Solid line has a slope 0.48 
in accordance with H14p. 

Figure [5] shows that the dependence of po (a; corr) is weh 
approximated by 

PoXcorr = Const ~ 0.48 (14) 

for different values of parameters of RQNLS and for 3 
different types of forcing. This gives another indication 
that the modulational instability determines the corre- 
lation length Xcorr through the amplitude of the back- 
ground fluctuations po in agreement with the equation 

m- 

Here, the first type is the standard forcing (O which, 
while not introducing any scale by itself (it pumps en- 
ergy into all Fourier modes fc), creates the scale Xcorr 
indirectly through the development of the modulational 
instability. Circles in Figure [5] correspond to that type 
of forcing. 

The second type is a particular example of the general 
type of forcing ^ . We define it by introducing a cutoff 
wave number kcutof / for the amplification; 

&fc = &o 7^ 0, |fc| < kcutof f\ h = 0, |fc| > kcutof f- 

(15) 
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FIG. 6: Time dependence of the number of particles N in 
the case of "fc-limited" deterministic forcing (f6)). (|15p with 
kcutof f ~ IOtt. Parameters are a = c = l,6o = 10*, e = 
2 • 10"^ 
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FIG. 7: Time dependence of the number of particles A'' in 
the case of random forcing. Parameters are a = c = l,e = 
10-^/c = 0.02,69 = 12500. 



In Figures 15161 for this second type of forcing we use (|6]) 
and (|15p with kcutof f = IOtt which means that only 11 
modes are amplified. The resulting N{t) dependence in 
Figure IHl is similar to Figure [5] In simulation of Figure 

El kjYiax — — 20 ^ kmaxi and thus, kcutof f ~ kmax- 

Squares in Figure [S] correspond to this type of forcing. 

The third type is a random additive forcing ([S]). Fig- 
ure [7] shows N{t) dependence from the simulations with 
a random forcing and the result is similar to Figures [2] 
and [HI Diamonds in Figure [5] correspond to this type of 
forcing. 

We conclude that the turbulence in RQNLS is inde- 
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pendent of the type of forcing used (up to the normahza- 
tion). The own scale of the forcing (if exist) is not im- 
portant and the correlation length is determined by the 
modulational instability scale. Qualitative picture is the 
following: the interplay between forcing and dissipation 
in RQNLS defines the amplitude po and the correlation 
length Xcorr to be related by ((TH). 

We also note that the principle difference between 
RQNLS turbulence and Navier-Stokes turbulence is the 
absence for RQNLS of well-defined inertial interval at 
which both forcing and dissipation are not important. 
For RQNLS dissipation and forcing generally act at all 
scales. Thus study of the spectrum of (IV'feP) could give 
much less information compare with e.g., weak turbu- 
lence [21*1. Instead below we focus on the study of the 
tail of PDF for large fluctuations. 



III. COLLAPSE AND ITS REGULARIZATION 
IN RQNLS 

Fluctuations of background result in multiple forma- 
tions of collapses in RQNLS as shown in Figure [TJ For 
e = 0, the collapsing solution of RQNLS has the following 
self-similar form 



rt ^t' 



r 



(16) 



where L — L{t) is the dynamically evolving spatial scale 
of the collapsing solution at a given moment of time, p 
and r are sometimes called by blow up variables. V{p, r) 
is well approximated for small i ^ 1 by the ground state 
soliton solution i?o(p) : ^(Pi ''') — Ro{p)- Here Rq{x) 
is the positive-definite solution of the equation — i?o + 
(9^i?o + -Rq — 0: which follows from RQNLS assuming 
that '0 = e^^Ro{x) and e = 0. The explicit expression for 
Rq{x) is given by 



Roix) 



31/4 



(cosh 2x) 1/2 



(17) 



The Hamiltonian Q (/? = 1 according to ([5])) vanishes 
at the groimd state soliton solution ([T7| . The number of 
particles (fTT|) at (fT7| is given by 



(18) 



Nc determines the boundary between collapsing and non- 
collapsing solutions: collapse is impossible for N < Nc- 
Note that the ground state soliton solution for NLS ([T]) 
with l3 — in dimension two (2D) is often referred to as 
the Townes solution and it plays a similar role to (jl7p in 
the collapse of 2D NLS. As L{t) decreases with t — > to, 
the spatial distribution of V{p) — ?> R{p) for p ^ I, where 
to is the collapse time (time at which a singularity devel- 
ops in the solution of RQNLS with e = 0). It implies 
that the number of particles in the collapsing region, 



NcoUpase Nc a,s t to, i.e., the collapse of QNLS 
is strong one as opposed to the weak collapse [ij, HO] 
for which the number of particles in the collapsing region 
vanishes as t — > io- 

The leading order behavior of L{t) for t to can be 
estimated from the scaling analysis as L{t) ^ {to — ty^'^ 
but the criticality of ID QNLS results in the following 
log-log modification of that scaling [23, [3ll - l33| 



2n{to - t) 
lnln[l/(to-t)] 



1/2 



(19) 



RQNLS with e > does not allow singular col- 
lapses. Instead, for < e ^ 1, the collapse amplitude 
IV'lmax(t) = max |'(/'(x, t)| goes (as a function of time) 

through a maximum \^p\maxmax = max \tp\max{t) at some 

t = tmax and decays after that as shown in Figure |H^. 
The spatial form of the collapsing solution is still well 
approximated by and p7)) before and shortly after 
t = tmax as seen in Figures [9K,b and c. It is seen in 
Figure that with the growth of \tp\max{t) the normal- 
ized shape of solution approaches to (IT7)) . But with the 
decay of liplmaxit) for t > tmax, the normalized shape of 
solution departs from (1171) and produces growing oscillat- 
ing tails as seen in Figures [Hb and O;. In the vicinity of 
t = imax the forcing on the rhs of RQNLS (O can be ne- 
glected. The resulting equation can be written in rescaled 

units tlV'liaxmaa;. AMliaxmaxi and 1p / \ll^\maxmax tO have 

exactly the same form, i.e., RQNLS without forcing is 
invariant with respect to these scaling transformations. 
As shown in FigureHb vs. Figure[5^, \^\max{t) rescaled 
in these units exhibits a universal behavior: all curves 
collapse on a single curve in the neighborhood of large 
collapses. That universality is independent of the com- 
plicated structure of optical turbulence. We conclude 
up to rescaling all collapse events in RQNLS are iden- 
tical which is qualitatively similar to the universality of 
collapse in QNLS. This universality is a characteristic 
feature of the dissipative terms in RQNLS ([S]). 

A function 



-L- 



dL 
Itt' 



(20) 



changes slowly with t compared to i at t ^ tmax except 
very small neighborhood of < = tmax as shown in Figure 
[TOl The dependence of ^{t) is shown in Figure [TOl in the 
rescale time unit {t—tmax)\i^\maxmax similar to Figure[5)D. 
We determine L{t) numerically from individual collapse 
events with \'4'\maxmax > 30 using and p?)) as L{t) = 
3^^'^ /\i^\max{ty ■ Thus, the universality of each collapse 
event is seen for 7(i) also. 

We conclude in this Section that the collapse events in 
RQNLS turbulence are all universal ones after a proper 
rescaling. 



7 



(a) 



200 r 



150 



_E 100 




-°10 



t-t 



X 10 




500 



FIG. 8: (Color online) (a) Dependence of li/ilmai on (t — tmax)- 
75 individual collapse events are collected for ji/'lmaa; > 15 
with parameters a — l,b = 10*, c = 1, e = 2 • 10~^. (b) 
The same dependence as in (a) but in rescaled units {t — 
tmax)\i>\maxmax and \tp\max / \ii\maxmax ■ The time intervals 



over which Itplmax /\^\7naxmax are plotted correspond to the 
time intervals of (a) (i.e., temporal interval is fixed in non- 
rescaled units). As a result, collapses with relatively small 
\'4}\maxmax extcuds over small intervals of {t — tmax)\ii\ 
in (b). 
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IV. PDF OF Iv^l 
A. PDFs from simulations 

Once the amplitude of a collapse event |'0|ma2; reaches 
the maximum \ip\maxmax, it starts decreasing and sub- 
sequently, the collapse event decays into outgoing waves 
as seen in Figure [SJa and[3J:. Decaying waves are almost 
liner ones as can be seen from the time dependence of the 



FIG. 9: (Color online) Scaled profiles of numerical solution 
|i/'(a;)| in the vicinity oi t — tmax at which |'i/>|mai reaches the 
maximum li^lmaxmax- Parameters are e = 2 • 10~^,a = c = 
l,b — 10*. (a) Scaled numerical solution |'i/'(a;)| at t = tmax — 
(4 • 10"^) (dotted red), t = t^ax - 10"* (dashed blue), and 
t = tmax (dashed-dotted green). Solid black line shows the 
normalized ground state soliton R{p) = 3~^'^*-Ro(p) (flT)) . (b) 
Scaled numerical solution |V'(3;)I at i = tmax — 10~* (dotted 
red), t = tmax (dashed blue), t = tmax + 10~* (dashed-dotted 
black), t — tmax + (7 • 10~*) (solid green). Circles show the 
normalized ground state soliton as in (a), (c) Logarithmic 
scale of Figure (b). 
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FIG. 10: (Color online) Dependence of j{t) — —LtL on (t — 
tmax)\4>\maxmax- The Same individual collapse events as in 
Figure |8] are collected for |V'|maa^ > 30 in non-rescaled units. 



ratio of the kinetic energy K = J \dxip\'^dx and the po- 
tential energy P = — ^ J \ip\^dx in Figure [TTk. Note that 
the Hamiltonian ([3]) is a sum of K and P and during the 
growth of the amphtude |'0|maa; of collapse event, these 
terms nearly cancel each other if the integrals are calcu- 
lated only inside the collapsing region p < 1. While at 
the decay phase of the collapse event, the amplitude that 
balances is strongly violated in favor of K meaning that 
solution is becoming nearly linear one. Figure [TTb shows 
zoom-in at a single collapse event. It is seen that while 
amplitude of collapse event decays, the kinetic energy 
from outgoing waves dominates over potential energy. 
Superposition of many of these almost linear waves from 
multiple collapse event forms a nearly random Gaussian 
field by the central limit theorem 34] . That random field 
seeds new collapse events. 

Figure [T2k shows PDF Vr{hr) for the real part of the 
amplitude V' to have a value hr- PDF for the imaginary 
part of tfj has the same form. Solid line is the match to 
the Gaussian distribution (with the same variance as for 
Vrihr)) which is almost indistinguishable from Vr{hr). 
PDF Vr{hr) is determined from simulations as 



J 5{Re{iJj{x,t)) - hr)dxdt 
J dxdt 



(21) 



Here the integrals are taken over all values of x and all 
values of t after the turbulence has reached the statisti- 
cally steady state. We assume ergodicity of turbulence, 
i.e., that averaging over space and time is equivalent to 
the averaging over ensemble of initial conditions (or the 
stochastic realizations of random forcing for the random 
forcing case). 

In a similar way, V{h) for the amplitude |^| to have a 
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FIG. 11: (Color online) (a) Time dependence of y^, where 
K is the kinetic energy and P is the potential energy defined 
in the text. Parameters are a= l,h = 10**, c = 1, e = 2 • 10"'^. 
Thick solid curve (blue) corresponds to Thin solid curve 
(green) is time dependence of maximal spatial value of 
Scale for ^ is on the left vertical axis and scale for max \ ip\ 
is on the right vertical axis, (b) Zoom-in at the same curves 
as in (a) for a single collapse event. 



value h is determined from simulations as 



V{h) 



J SiltPix^t)] - h)dxdt 
J dxdt 



(22) 



and, as shown in Figure [T^b, it is again very close to the 
Gaussian distribution. Large fiuctuations of l^/^l are how- 
ever, quite different from the Gaussian distribution and 
have power-like tails as shown in Figure [T2h at log-log 
scale. From comparison of Figures [T^ . b and c we con- 
clude that the fit to the Gaussian distribution works very 
well for IV'I ^ 8 which can be interpreted as the super- 
position of numerous almost linear waves. For \^\ <; 10 
the PDF has a power law-like dependence which can be 
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FIG. 12: (Color online) PDFs for simulations with a — l,c = 
1,& = W.e = 2-10-^ (a) PDFVr{hr)iorhr = Re{ip). Solid 
blue line represents numerical result and dotted blue line is 
for Gaussian distribution (27r)"^''^/io ^e"''^/'^''"' plotted for 
comparison. Here the variance /iq = 5.85 is obtained from 
the simulation, (b) PDF of h = I'i/;! at linear scale. 

Solid blue line shows numerical result and dotted blue line is 
Gaussian distribution h^^he'*^ /(z^o)^ where /iq is the same 
variance as in (a), (c) PDF V{h) oi h = \%p\ &t log-log scale. 
Solid blue line shows numerical result and dotted blue line is 
the Gaussian distribution as in (b). Thick solid red line shows 
law for comparison. 




FIG. 13: (Color online) (a) PDF V{h) of ft = for "fc- 
limited" deterministic forcing ((6]), (|15|) with kcutoff = IOtt. 
Parameters are a = c = l,&o = 10"*, e = 2 ■ 10"^. (b) PDF 
V{h) of h = lipl for random forcing with parameters a = c = 
1, e = 10~^ = 0.02, bg = 12500. 



roughly estimated a.s ^ h ^ and which indicates inter- 
mittency of the optical turbulence (25| . 

Figure [T3ji shows the PDF r{h) for RQNLS with the 
"fc-limited" deterministic forcing ([5]), (ITSt and Figure [T^ 
shows the PDF V{h) for RQNLS with the random forcing 
It is seen that both of these cases have the same type 
of power-like tails as in Figure [T2b. Thus, the power-like 
tails are universal for RQNLS for e ^ 1 and independent 
of the type of forcing indicating the universal turbulent 
behavior. 

In order to characterize the relevance of linear and non- 
linear dissipation in PDF V{h), we consider the different 
parameter values for a and c while fixing e and b. Fig- 
ure 114b . shows V{h) for fixed nonlinear dissipation c — 1 
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and different values of the linear dissipation coefficient 
a. As a decreases, the slope of PDF tail becomes steeper 
deviating from power law. For a < 0.1, we observe 
that the shape of V{h) tail is very close to the case of 
a = 0.1. Figure [Mb shows for fixed linear dissipa- 

tion a = 1 and variable c. The PDF tail again deviates 
from h^^ power law as c decreases. The change of tail 
is more significant at large amplitudes h. As c becomes 
smaller, high amplitude collapses occur more often, e.g., 
for c ~ 0.05, we frequently observe collapses with ampli- 
tude > 500. In contrast, if we change e only, then the 
tails of Vih) behave similarly for small e ^ 3-10"^ except 
very large values of h as shown in Figure fT4b. 

We also study the sensitivity of 'P{h) to the change of 
amplification amplitude b. Figure [T51 shows that a height 
of the maximum of Vih) decreases with the increase of 
b while the position of the maximum h = Hq shifts to 
the right indicating the increase of the average ampli- 
tude (IV'P) (see also Figure[5]). To stress this feature we 
focus in Figure [15] on a smaller domain in h and com- 
pare with previous figures, while for larger h the tails of 
V{h) behave similar for different b up to the normaliza- 
tion constant. 



B. Power law-like tails of PDF 



We now show that the power-like tail of 'P{h) results 
from the near-singular collapsing events. This approach 
dates back to the idea of describing strong turbulence in 
the Navier-Stokes equations through singularities of the 
Euler equations [2^ . Unfortunately, this hydrodynamic 
problem remains unsolved. The forced Burgers equation 
remains as the only example of an analytical description 
of strong turbulence in which the tail of the PDF for neg- 
ative gradients follows a well established (—7/2) power 
law [23|, dominated by the dynamics of near-singular 
shocks. Another example of the analytical description of 
the intermittency is a randomly advected passive scalar 
which is the example of the turbulent transport described 
by linear equations [s^ . 

As a first step we calculate the contribution to the PDF 
from individual collapse events. As shown in Figure |S1 
the filament amplitude l-iAlmaa: rapidly decays after reach- 
ing {iplmaxmax at < = tmax- Thus, wc ueglcct the Contri- 
bution to V{h) from t ^ tmax in calculating the con- 
tribution of the individual filament to V{h). We define 
the conditional probability V {h\hmax) for contribution 
to PDF from the collapse event with \ip\maxmax = hmax 



(a) 




h 



FIG. 14: (Color online) PDF T(h) oih^ (a) Parameters 
are e = 2 • 10"^, c = 1,& = 10*, a = 0.05 (dotted), a = 
0.1 (dashed), a — 0.5 (dashed-dotted), a = 1 (solid), (b) 
Parameters are e — 2 ■ 10~^, a — l,b — 10*, c — 0.1 (dashed), 
c = 0.5 (dashed-dotted), c = 1 (solid), (c) Parameters are 
a = l,c = 1,6 = 10*, e = 10"^ (solid), e = 2 ■ 10"^ (dashed- 
dotted), 6 = 5- 10~^ (dashed). Thick solid line shows 
power law in all Figures. 
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FIG. 15: (Color online) Zoom-in of PDF V{h) on the range 
of small h = \ip\- Parameters are e = 2 ■ 10"'^, a = c — 1, 
& = 2 ■ 10^ (solid), b = 10* (dashed), 6 = 5 • 10* (dashed- 
dotted) . 

and use ([HI) ,([I71), (Uni and dH]) as follows 

= Const h <d {hmax - h) , 

where h^aax = M^) I L{tmaxY'^ = ^-^^ ^ I L{tmaxf , and 
Q{x) is the Heaviside step function. Here, we have 
changed the integration variable from t to L and approx- 
imated j{t) under the integral by its average value (7) 
as "f{t) ~ (7) ~ 0.5. This approximation is valid for 
t ^ tmax outside the neighborhood oft — t^ax as seen in 
Figure Uni 

As a second step we calculate 7^(/i) by integration over 
all values of h^ax using equation (|23p as follows 

^(^) j dhYy^Q^x^ih\hmax^^rnax{hraax^ 

— Const h dhjjiQx^i^max ^^T^maxi^max^ 

= Const h~^Hmax{h), (24) 

where l^raaxij^max^ the PDF for h^y^dx — l^lmaa^maa: 

and Hmaxih) = Vrnax{hrnax)dhrnax IS thc Cumulative 

probability for \tl^\maxrnax > h. 

Figures 116b and 116b show Hmaxihmax) for different 
values of parameters. Each curve is calculated after the 
system reaches the statistical steady using more than 



10'^ collapse events with hmax > 10. We verified that 
the increase of the number of collapse event (i.e., in- 
crease of the total simulation time) does not change 
these curves in any significant way. Figure 116b , shows 
that power-like dependence h~^^x for Hmax{Knax) if ex- 
ist at all, quickly disappear with the increase of e. Fig- 
ure [TBJd shows that Hmaxihmax) shows that power- like 
dependence ft'^^ shifts to larger values hmax with the 
decrease of a. Generally, we see from Figures [T6k and 
[T5b that for a wide range of parameters, including the 
decrease of e, the dependence of Hmaxihmax) cannot 
be approximated as oc (hmax)~^ ■ The assumption of 

Hmax {hmax max 

) ^ comes from and the rough 
estimate that V{h) oc h^^ as in thick solid red line of 
Figure [T^. We conclude that this conjecture, first made 
in ,27.] . appears to be incorrect and is only a very 
crude approximation for Vih). 

We now verify that the equation (1^^ is correct one, 
which justifies the assumptions used in (|23p and (|24l) . 
Figures [TTk. [TTb andfTTb compare Vih) from simulations 
(solid blue lines) with the prediction of the equation ([M)) 
(red circles), where Hmaxihmax) is obtained numerically 
and shown in Figure [1^1 Also Figure [1^] shows that for 
h ^ 20 the equation (|24l) appears to be much better fit 
of Vih) compare with h^^ power law. The normalization 
constant in the equation ([M)) was chosen to fit Vih) at 
large h. The deviation of (HH) from Vih) for /i < 20 is 
perhaps due to the decaying of the large amplitude col- 
lapse events into large amplitude waves (as e.g. in Figure 
[5J;). This explanation is consistent with the increase of 
that deviation for smaller values of a as in Figures [T7b 
because the decrease of a causes these large amplitudes 
almost linear waves to live longer before dissipating, thus 
giving a bigger contribution into Vih). We also conclude 
from Figures [TTk . [TTb and [TTb that because such devia- 
tion is insignificant for h > 20, the contribution of large 
amplitude waves is not important in that range of h. 
The good agreement between Vih) and the equation ([M]) 
justifies the assumptions used in derivation of the equa- 
tion as well as it shows that the intermittency of 
optical turbulence of RQNLS (O is solely due to collapse 
dynamics. We also conclude that while h~^ power law 
appears to be an intermediate fit at best, the equation 
([M)) works for all values of parameters of RQNLS we 
tested. 



V. NUMERICAL SIMULATIONS OF ID RQNLS 

In this section we provide detailed description of the 
numerical methods used to produce the simulations de- 
scribed above. We conduct numerical simulation of 
RQNL ([5]) by employing a version of the fourth order 
split-step method [36| (outlined below) for deterministic 
forcing as well as the second order split-step method for 
random forcing. For efficient computation, we adaptively 
change the spatial grid size Ax during the time evolution 
in the spatial domain —0.5 < x < 0.5. Specifically, if the 
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FIG. 16: (Color online) Cumulative probability Hn-iax{hmax j 
for hmax = \ip\max soid deterministic forcing with c = l,b — 
10*. (a) Parameters are a = l,e = 5-10"^ (solid), e = 2-10"^ 
(dashed-dotted), e = 10"'^ (dashed), (b) Parameters are e — 
2 • 10"^. a = 0.05 (solid), a = 0.1 (dashed-dotted), a = 0.5 
(dashed), a = 1 (dotted). Thick solid red line shows h. 
power law for comparison in both figures. 



1 

max 



amplitudes of Fourier components of solution, \ijjk\, at 
high frequencies exceed lO"^*' max {tpkl, we reduce Ax by 

k 

adding more Fourier modes to the system. If the am- 
phtudes of high frequency modes are below the criteria, 
we increase Ax by removing some of the existing Fourier 
modes. The numerical time step. At, is also updated as 
Aa; changes, which follows the relation Atk^^^ — qoir. 
Here k^ax — tt/Ax is the maximum wavenumber de- 
termined by the discretization and we choose the con- 
stant factor Qo small enough to avoid numerical insta- 
bility. The numerical instability occurs if the change of 
phase A(j) = Aifc^^^ of the highest Fourier harmonics 




FIG. 17: (Color online) Solid curve represents PDF P(/i) 
for h = \ip\, and dashed line shows power law. Cir- 
cles correspond to the rhs of equation (|24l) . (a) Parame- 



ters are e 
ters are e 
e = 2 ■ 10" 



= c = l,b = 
10~^a = c = l,fe = 10". 
^a = 0.5, c= 1,6= 10*. 



2 ■ 10" 

-3 



10*. (b) Parame- 
(c) Parameters are 
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kmax from the linear term of RQNLS at one discrete time 
step At is above 7r/2 : A(f> > tt/2. In that case a mixing 
of the Fourier harmonics in the quintic nonhnear term in 
RQNLS can produce artificial (discretization-caused) res- 
onance from the condition AAcj) = 27r. It implies that to 
avoid the instability we have to choose qq < 1/2. In simu- 
lation we typically choose qo = 0.2 which allows to avoid 
instability as well as insures high accuracy in time step- 
ping. The initial condition of our numerical simulation 
is a random field with small amplitude whose maximum 
value is about 1.6. 



A. Fourth order Split-Step scheme 

We write equation ([5|) formally in the form of 



(25) 



where L is the linear and N is the nonlinear operators, 
defined by 



Lip = i{l - iae)\'^ip + e<p, 
Nip = i{l + ice)\ip\^il;. 



(26) 
(27) 



Separately both ([26)) and (|27)) can be solved very 
efficiently. Here, we approximate the exact solution 

ip{to + At) = e^*(^+^V(io) of (HSI over one segment 
from to to to + At by the following expression 



1p{to + At) = eC4AtAfgd3AtLgC3AtA'gd2AtLgC2AtJV 



(28) 



where 



ci = 



1 

2(2-21/3) I 
C3 = C2, 



C2 
C4 



2(2-21/3) ' 



di = 



1 

2-21/3 ■ 



_2i/3 
2-21/3 , 



di. 



(29) 



The split-step expression (1^51) . is of the fourth or- 
der accurate and it is a straightforward generalization 
of the forth order symplectic integration of [361 non- 
Hamiltonian systems. 

The linear part e'^*^ of the operator splitting can be ef- 
ficiently calculated by using the Fast Fourier Transform 
algorithm (FFT) for the deterministic forcing ([6]). For 
the stochastic forcing ([8]) the equation ([26l) has a form 
of the inhomogeneous linear differential equation. Ho- 
mogeneous part of that equation we solve again using 
FFT while the contribution of the inhomogeneous term 
is obtained by the numerical integration over t for each 
X. We use the trapezoidal rule for the integration which 
makes the scheme the second order in the random noise 
case. Therefore, instead of the fourth-order split-step al- 
gorithm , we use the standard second order split-step 
algorithm ^(to + Ai) = e(i/2)AtigAt7Vg(i/2)AtL^(^^)^ 
the random noise case. 



The nonlinear part e^*^ of the operator splitting can 
be solved exactly as follows. Denote the solution of the 
nonlinear part of RQNLS as ip^ . It means that V'^ (^o + 
At) = e^**i/;*(to). According to RQNLS, ip'^ satisfies 
the following equation: 



dt^^ = {-ce + t)\tP^\^iP^. 



It implies that 



dt 



N\2 



-2ce|V'^|6, 



(30) 



(31) 



and hence, we find 

1^-^(^)12=0, if |V*(io)P = 0, (32) 
\iP^{t)\^ = [4ce(f- to) -K|i^*(to)n-i otherwise. 

(33) 

Using equations ([32 ]) . (|33l) . we obtain the explicit solution 
of equation (|30p . 



^''{to + At) 



exp 



X V'^(to) 



4ce 



■\n[AceAt\^^{to)f + l] 

(34) 



In the case of random forcing (jS]), at each time step 
we independently generate the random variable ^(x) by 
the Ornstein-Uhlenbeck process ^3] as a function of x 
with zero mean and correlation length Ic = 0.02 that sat- 
isfies a stochastic differential equation, ^ = —l^^^{x) + 

a{At)-^/^^^^, where ct > and W{x) denotes the 
Wiener process. This Ornstein-Uhlenbeck process yields 
an exponential correlation function (|9]) with bg — 
Here the factor (Ai)^^/-^ multiplying dW/dx ensures 5- 
correlation of ^ in time for At — > 0. 



VI. CONCLUSION 

In this paper we studied the strong turbulence in ID 
RQNLS ([5]). In the statistical steady-state (the state of 
developed turbulence) the dynamical balance is achieved 
between forcing which pumps the number of particles 
N into the system and both linear and nonlinear dissi- 
pation. RQNLS has multiple collapse events randomly 
distributed in space and time. We found that in the 
state of developed turbulence the spatial and tempo- 
ral correlation functions have universal forms indepen- 
dent of the particular values of parameters of RQNLS as 
well as independent of the type of forcing. In particu- 
lar we considered the deterministic forcing with variable 
fc— dependence as well as random forcing. We found that 
the relation between the correlation length and the av- 
erage amplitude po has a universal form (jl4p determined 
by the modulational instability scale. 

PDF V{h) of amplitude fluctuations is well approxi- 
mated by the Gaussian distribution for h ^ 3po- In con- 
trast, 'P{h) for h ^ 4po has the strongly non-Gaussian tail 
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with power-like behavior characterizing intermittency of 
strong collapse-dominated turbulence. This tail is deter- 
mined by the equation (|24|) which includes the contri- 
bution (X /i"*" from the universal spatio-temporal form 
of the collapse events as well as the contribution from 
the cumulative probability Hmaxih), the probability of 
the maximum amplitude of collapse event that exceeds 
h. We show that Hmax{h) is not universal and depends 
on the parameters of RQNLS. For some range of param- 
eters Hjnaxih) can be roughly estimated as oc but it 
appears to be an intermediate asymptotic at best 

An important problem to be studied in future work 
is to determine the analytical form of H„iax {h) from the 
parameters of RQNLS. This is a challenging problem and 



will require calculation of the optimal fluctuations of the 
background which seeds new collapses. Moreover, it will 
be necessary to find the relation between the size of such 
optimal fluctuation and the maximum collapse event am- 
plitude \i^\maxinax- 
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